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Introduction 

Considerable work has been done in recent years 
on the development of contact algorithms for use in 
conjunction with finite element programs. These ac- 
tivities were motivated by a variety of practical ap- 
plications including investigation of wear resistance 
of machine elements in contact, determination of the 
load distribution in joints, and determination of the 
footprint area of tires. Recent work on this subject 
is reviewed by Kikuchi and Oden (1988); Noor and 
Tanner (1988); Simo, Wriggers, and Taylor (1985); 
Torstenfelt (1983); and Wriggers, Wagner, and Stein 
(1987). 

Contact algorithms first used in finite element ap- 
plications were based on semiempirical approaches 
such as the gap element approach (see, for exam- 
ple, Stadter and Weiss 1979). These approaches of- 
ten lacked the generality needed for handling com- 
plicated practical problems. Since contact conditions 
are formulated as inequality constraints, mathemat- 
ical programming approaches were naturally used to 
incorporate the contact constraints into the global 
finite element equations. The two most commonly 
used methods in this category are the Lagrange mul- 
tiplier technique and the penalty method (see, for ex- 
ample, Felippa 1978; Kalker 1986; and Kikuchi and 
Song 1981). The main drawbacks of the Lagrange 
multiplier technique are that it results in an increase 
in the number of algebraic equations to be solved 
and the presence of zero terms along the main diag- 
onal of the coefficient matrix of these equations. In 
the penalty method the contact conditions are only 
approximately satisfied and the resulting algebraic 
equations become ill-conditioned for large values of 
the penalty parameter. 

More recently, a perturbed Lagrangian formula- 
tion was proposed for the solution of contact prob- 
lems. The formulation includes the Lagrange multi- 
plier technique and the penalty method as limiting 
cases and alleviates some of their drawbacks. Al- 
though the perturbed Lagrangian formulation can 
be interpreted as a mixed method, it has been used 
only in conjunction with displacement finite ele- 
ment models (see Simo, Wriggers, and Taylor 1985; 
Stein, Wagner, and Wriggers 1986; and Wriggers, 
Wagner, and Stein 1987). Mixed formulations (with 
the fundamental unknowns consisting of the gener- 
alized displacements and internal forces) were shown 
to be considerably simpler for the treatment of highly 
nonlinear problems than displacement formulations. 
Therefore, the question arises as to whether the use 
of the perturbed Lagrangian formulation in conjunc- 
tion with mixed finite element models might result in 
a more effective computational strategy for contact 


problems. The present study focuses on this ques- 
tion. Specifically, the objectives of this paper are 
(1) to present simple mixed finite element models and 
a computational procedure b8ised on the perturbed 
Lagrangian formulation for the solution of contact 
problems and (2) to demonstrate the effectiveness of 
the models and the proposed procedure by means of 
numerical examples. 

To sharpen the focus of the study, only frictionless 
contact of axisymmetric shells and curved beams is 
considered. The effects of large rotations and trans- 
verse shear deformation are accounted for, and the 
material of the structures is assumed to be linearly 
elastic and isotropic. 


Notation 

A 

E 

[F] 


f{Z,p) 

{G{Z)} 


9o 

{9o} 


cross-sectional area of the beam 

Young’s modulus of the material 

flexibility matrix for an individual 
element 

shear modulus of the material 

vectors of nonlinear terms in the 
element equations 

vector defined in equations (6) 

vector of nonlinear contributions to 
the global equations 

current gap (measured in the 
direction of the normal to the 
contact surface) 

initial gap 

vector of initial gaps for the contact 
element 


{H} 


vector of internal force parameters 


h shell thickness 

/ moment of inertia of the beam cross 

section 



Ms.Mq 


N 


global linear matrix of the structure 

moment resultant turning about the 
normal to the middle surface of the 
shell (see fig. 1) 

meridional and circumferential 
(hoop) bending stress resultants 
(see fig. 1) 

shape functions used for approx- 
imating the generalized displace- 
ments and the Lagrange multipliers 
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Na^ No meridional and circumferential 

(hoop) stress resultants (see fig. 1) 

P concentrated load on beam (or ring 

load on shell) 

{P} normalized external load vector 

{P} global vector of normalized external 

loads and initial gaps 

Pn nodal force normal to the contact 

surface 

p load parameter 

Po intensity of pressure load (internal 

pressures have positive sign) 

[Q], [i?] elemental matrices associated with 

the contact condition and the 
regularization term in the functional 

Qa transverse shear stress resultant (see 

fig- 1) 

r radius of hemispherical shell (see 

fig- 1) 

Vo radius of curvature of circular ring 

and hemispherical shell 

[S'] strain-displacement matrix for an 

individual element 

s tangential coordinate along the 

undeformed centerline of the beam 
or meridional coordinate of the 
undeformed shell (see fig. 1) 

Se values of s at the initial and the last 

point of the contact region 

Tfi intensity of contact pressure (acting 

normal to the contact surface) 

U strain energy of the structure 

U strain energy density (strain energy 

per unit length for beam or per unit 
area for shell) 

u, w displacement components in the x\ 

and X 3 coordinate directions (see 

fig- 1) 

{X} vector of nodal displacements 

^3 Cartesian coordinates (for a shell 

X 3 coincides with the axis of revolu- 
tion; see fig. 1) 

{Z} global response vector 


preselected constants used in the 
constraint equation (see eqs. (7) 
and (8)) 

7 o transverse shear strain 

e penalty parameter 

Ss^ja virtual extensional and transverse 

shear strains 

Sao extensional strain (in the meridional 

direction of the shell or along the 
length of a beam) 

sq extensional strain in the circumfer- 

ential direction for a shell 

C generalized arc-length parameter 

defined in equations (7) and (8) 

Kji bending strain associated with the 

moment turning about the normal 
to the shell Mn 

Ks bending strain for a beam (or 

meridional bending strain for a 
shell) 

K 0 circumferential (hoop) bending 

strain for a shell 

A Lagrange multiplier, representing 

the intensity of the contact pressure 
(acting normal to the contact 
surface) 

{A} vector of nodal values of the La- 

grange multiplier 

u Poisson’s ratio of the material 

n functional 

(f) rotation component of the beam or 

the middle surface of the shell 

(f> = 0 + 

00 angle which a typical cross section 

of the undeformed beam makes with 
the a; 3 -axis (or angle between the 
axis of revolution and the normal to 
the shell middle surface, see fig. 1) 

element domain 
fie contact surface 

da = djds 

Subscripts: 

a at point a of circular ring in 

figure 11 
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b at point b of hemispherical shell in 

figure 3 

conv, converged 

Superscript t denotes transposition. 

Finite elements are designated M3-2 and M4-3 as 
shown in table 1. 

Mathematical Formulation 

The analytical formulation is based on a form of 
Reissner’s large-rotation theory with the effects of 
transverse sheetr deformation included. A mixed for- 
mulation is used in which the fundamental unknowns 
consist of the three generalized displacements and 
the internal forces (five for shells of revolution and 
three for curved beams). The sign convention for the 
generalized displacements and the internal forces is 
given in figure 1. The fundamental equations of the 
axisymmetric shell theory and the curved beam the- 
ory used herein are given by Reissner (1972a, 1972b) 
and are summarized in appendix A. 

Contact Condition 

Figure 2 shows the characteristics of the friction- 
less contact of a structure pressed against a rigid 
plate: ilc refers to the contact region; is the initial 
gap between the structure and the plate; g is the cur- 
rent gap (both g^ and g are measured in the direction 
of the normal to flc); and Tn is the normal traction 
on He- The contact condition can be expressed by 
the following inequalities, which must be satisfied at 
each point on the contact surface flc: 


contact pressures and a regularization term which 
is quadratic in the Lagrange multipliers. For a de- 
tailed discussion of the perturbed and augmented 
Lagrangian formulations, see Simo, Wriggers, and 
Taylor (1985); Stein, Wagner, and Wriggers (1986); 
and Wriggers, Wagner, and Stein (1987). 

The modified functional has the following form: 

n = Ujjji + ~ ^ (^) ] ^ 

where Hhr is the funcUonal of the Hellinger- Reissner 
variational principle; A is the Lagrange multiplier; 
and e is a penalty parameter associated with the reg- 
ularization term. The explicit form of for planar 
deformation of curved beams is given by Noor, Pe- 
ters, and Andersen (1984); and for axisymmetric de- 
formation of shells of revolution, by Noor, Andersen, 
and Tanner (1984), Note that the addition of the reg- 
ularization term amounts to approximating the rigid 
plate by continuously distributed springs with stiff- 
ness of e, for sufliciently large e. As e approaches 
zero, the continuous springs become the rigid plate. 

The shape frmetions used in approximating the 
generalized displacements and the Lagrange multipli- 
ers are selected to be the same and differ from those 
used in approximating the internal forces. More- 
over, because of the nature of the functional n in 
equation (2), the continuity of both the internal 
forces and the Lagrange multipliers is not imposed at 
interelement boimdaries. 

The finite element equations for each individual 
element can be cast in the following compact form: 


9>0 Tn<0 Tng = 0 ( 1 ) 

The first inequality in equations (1) represents the 
kinematic condition of no penetration of the rigid 
plate (zero gap for the contact points); the second 
inequality is the static condition of compressive (or 
zero) normal tractions; and the third condition is 
that of zero work done by the contact stresses (i.e., 
the contact stresses exist only at the points where 
the structure is in contact with the rigid plate). The 
inequalities p > 0 and Tn> 0 are henceforth referred 
to as the “inactive contact conditions.” 

Governing Finite Element Equations 

The discrete equations governing the response of 
the structure are obtained by applying a modified 
form of the two-field Hellinger-Reissner mixed vari- 
ational principle. The modification consists of aug- 
menting the functional of that principle by two terms: 
the Lagrange multiplier associated with the nodal 




r G{x) ] 

(«) r . 'I 

+ i 

M{H,X) 

► - IpP 


1 • J 

( 9o . 


( 3 ) 


where {-ff}, {A"}, and {A} are the vectors of the in- 
ternal force parameters, nodal values of the general- 
ized displacements, and nodal values of the Lagrange 
multipliers; [F] is the matrbe of linear flexibility co- 
efficients; [5] is the strain-displacement matrix; [Q] 
and [R] are matrices associated with the contact con- 
dition and the regularization term in the functional 
(see appendix B); {G(A)} and {M{H,X)} are vec- 
tors of nonlinear terms; {go} is the vector of initial 
gaps in the contact region; a dot refers to a zero 
submatrix or a zero subvector; superscript (e) refers 
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to individual elements; {P} is the normalized exter- 
nal load vector; and p is a load parameter. As the 
load is incremented, only the value of the load pa- 
rameter p changes and the normalized vector {P} 
is constant. The formulas for the elemental arrays 
[P],[5],{G(A)},{M(P,X)}, and {P} are given by 
Noor, Peters, and Andersen (1984) for curved beams 
and by Noor, Andersen, and Tanner (1984) for axi- 
symmetric shells. The formulas for the elemental ar- 
rays [Q] and [R] are given in appendix B. 

Note that the size of the coefficient matrices 
[P], [Q], and {go} varies with the number of ac- 
tive contact conditions. The difficulty associated 
with an equation system whose size varies during 
the solution process can be alleviated by allow- 
ing the Lagrange multipliers to be discontinuous at 
interelement boundaries and eliminating them on the 
element level. This was done in the present study. If 
the internal force parameters and the Lagrange mul- 
tiplier parameters are eliminated from equations (3), 
one obtains the following equations in the nodal dis- 
placements {X}: 

[[S)‘[F|-*(S) - £|Q1M-‘I<3I']"’ 

-p{P}^^^=0 (4) 

where 

+ {M{H,X)} (5) 

and the vector {H} in [M{H,X)} is replaced by its 
expression in terms of {X}. 

Solution of Nonlinear Algebraic Equations 

The discrete equations governing the response of 
the entire structure are obtained by assembling the 
elemental contributions in equations (3) or (4) and 
can be written in the following compact form: 

{f{Z,p)} = [k]{Z} + {G{Z)}-p{P} = 0 (6) 

where [K] is the global linear matrix of the structure; 
[G{Z)} is the vector of nonlinear contributions; {P} 
is the global vector of normalized external loads 
and initial gaps; and {Z} is the global response 
vector of the structure (obtained by assembling the 
contributions from the subvectors {i/},{A}, and 

{A}). 

The nonlinear algebraic equations (eqs. (6)) are 
solved and the contact region and the contact 


pressures are determined by using an incremental- 
iterative technique (i.e., a predictor-corrector contin- 
uation method) in which the response vector {Z}, 
corresponding to a particular value of the load pa- 
rameter p, is used to calculate a suitable approxi- 
mation (predictor) for {Z} at a different value of p. 
This approximation is then chosen as an initial es- 
timate for {Z} in a corrective iterative scheme such 
as the Newton-Raphson technique. In each Newton- 
Raphson iteration the contact conditions are checked 
and updated. 

Generalized Arc Length and Associated 
Constraint Equation 

A number of numerical experiments have demon- 
strated that a more effective computational algo- 
rithm results if {X}, {i/}, {A}, and p are considered 
functions of a nondecreasing parameter ( which is de- 
fined by the following constraint equations (see, for 
example, Ramm 1981): 

(dHYidH] (dxY^dX} 

where a^f, ax, and (3 are preselected constants 
(0 < OL}{^ax,OLx,0 < 1), and superscript t denotes 
transposition. If ajj = ax = a\ = f3 = then 
C represents the arc length in the space spanned by 
{i/}, {X}, {A}, and p. Therefore, ^ in equations (7) 
represents a generalized arc length in the solution 
space. 

In the present study a^, c^x^ were all 

taken to be equal, that is, ajj = ax — a\ = a, so 
that equations (7) can be written in the form: 

To retain the symmetry and bandedness of the re- 
sulting finite element equations of the structure, the 
constraint equations (eqs. (7) or (8)) are not solved 
simultaneously with equations (6); rather, the two- 
step procedure outlined by Ramm (1981) is used (see 
appendix C). The use of C as the parameter control- 
ling the progress of computation along the solution 
path, in addition to circumventing the problems as- 
sociated with the singularity of the system matrix 
at limit points, has other computational advantages 
that are discussed in the section on numerical studies. 
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Computational Procedure and 
Determination of Contact Pressures 

The computational procedure used in determin- 
ing the contact region and the contact pressures is 
outlined in appendix C. The nonlinearities due to 
the large rotations and the contact condition are 
combined into a single iteration loop. Wriggers and 
Nour-Omid (1984) advocated the use of a two-level 
(nested) iteration scheme. The inner iteration loop 
accounts for the contact conditions associated with 
the contact pressures, and the outer iteration loop 
uses Newton-Raphson iteration. Numerical experi- 
ments have demonstrated that for frictionless con- 
tact problems the two-level iterative scheme requires 
more iterations than the single-level scheme. This is 
discussed further in the section on numerical studies. 

The solution of the governing discrete equations 
of the entire structure generates the nodal displace- 
ments, the internal force parameters, and the values 
of the Lagrange multipliers at the contact nodes. For 
each individual element in contact, the intensity of 
the contact pressure at a node Tn is equal to the 
value of the Lagrange multiplier A at the same node. 
The contact pressures are also related to the nodal 
forces normal to the contact as follows: 

dQ tI (9) 

" Jnie) ” 

where are the shape functions used in approxi- 
mating the Lagrange multiplier and the generalized 
displacements and is the domain of the con- 
tact element. The range of both i and j in equa- 
tions (9) is from 1 to the number of displacement 
nodes. Torstenfelt (1983) discusses other approaches 
for determining the contact pressures. 

Comments on the Mixed Models, 

Perturbed Lagrangian Formulation, 
and Computational Procedure 

The following comments regarding the mixed 
models, the perturbed Lagrangian formulation and 
the computational procedure used herein are in 
order: 

1. The nonlinear terms in the finite element 
equations of the mixed model (eqs. (3)) have a 
considerably simpler form than those of the cor- 
responding displacement model (eqs. (5)). This 
is particularly true for large rotation problems 
(for a more detailed comparison between mixed 
and displacement finite element models, see Noor, 
Peters, and Andersen 1984; and Noor, Andersen, 
and Tanner 1984). 


2. Equations (3) include both those of the Lagrange 
multiplier approach and the penalty method as 
special cases, as follows: 

a. By letting the penalty parameter e go to in- 
finity, equations (3) reduce to those of the La- 
grange multiplier approach. 

b. By eliminating the Lagrange multiplier terms 
from equations (3), the resulting equations are 
identical to those of the penalty method. 

3. The perturbed Lagrangian formulation alleviates 
two of the drawbacks of the Lagrange multiplier 
approach and the penalty method, namely, 

a. The regularization term in the functional re- 
sults in replacing one of the zero diago- 
nal blocks in the discrete equations of the 
Lagrange multiplier approach by the diagonal 
matrix [R]/e in equations (3). 

b. The contact condition is satisfied exactly by 
transforming the constrained problem to an 
unconstrained one through the introdwtion 
of Lagrange multipliers (the term jQ^Xgdil 
in eq. (2)) rather than approximately as in 
the penalty method. However, the pres- 
ence of the_ regularization term (the term 
—Jq (l/2e)(A)^ dQ in eq. (2)) results in re- 
placing the contact condition by the perturbed 
condition, 

i[/?]{A} + [Q]‘{X}-{ffo} = 0 (10) 

4. An important consideration in the perturbed 
Lagrangian formulation (and in any penalty for- 
mulation) is the proper selection of the penalty 
parameter e. With the foregoing mixed mod- 
els the penalty parameter can be chosen inde- 
pendently of the element size, without adversely 
affecting the performance of the model. The 
accuracy of the solution increases with increas- 
ing value of the penalty parameter. However, 
for very large values of e, the equations become 
ill-conditioned and thus round-off errors increase 
(see Felippa 1978). In the present study the 
penalty parameter e was related to the exten- 
sional stiffness of the structure. 

5. The elemental arrays [F], [5], {G(A)},{M(i/, A)}, 
and {P} are evaluated numerically using a Gauss- 
Legendre formula. The arrays [Q], [P], and {go} 
are evaluated using a Newton-Cotes formula. In 
both cases the number of quadrature points used 
is the same as the number of displacement nodes 
in the element. This results in underintegrating 
the arrays [Q] and [P] and avoids the oscillatory 
behavior of the contact pressures that has been 
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observed when these arrays are fully integrated. 
Note that the use of Newton-Cotes formula al- 
lows the contact pressures to be evaluated at the 
displacement nodes. 

Numerical Studies 

To test and evaluate the performance of the mixed 
models and the foregoing computational procedure, 
several frictionless contact problems have been solved 
by using this procedure. For each problem, compari- 
son was made with converged solutions and with pre- 
viously published numerical and analytic solutions 
(whenever available). Herein the results of two typ- 
ical frictionless contact problems are discussed. The 
two problems are (1) a hemispherical shell in con- 
tact with a rigid plate (fig. 3) and (2) a circular ring 
pressed against a rigid plate (fig. 11). Updike and 
Kalnins (1972) give an analytic solution for the first 
problem. Finite element solutions using displace- 
ment models are presented by Stein, Wagner, and 
Wriggers (1986) for the hemispherical shell and by 
Simo et al. (1984) for the circular ring. In the cited 
references only load-deflection curves and distribu- 
tion of contact stresses were presented. The material 
for both structures is assumed to be linearly elastic 
and isotropic, and the effects of both large rotations 
and transverse shear deformation are accounted for 
(see appendix A). For the hemispherical shell, only 
axisymmetric deformations are considered and there- 
fore, only one meridian was modeled. For the circular 
ring only planar deformation is considered, and be- 
cause of symmetry, only half the ring was modeled. 

In each case the structure was analyzed using the 
mixed models developed herein. Lagrangian inter- 
polation functions were used for approximating each 
of the generalized displacements, internal forces, and 
Lagrange multipliers. For convenience, the values of 
a and /? in the generalized arc-length constraint equa- 
tions (eqs. (6) and (7)) were selected to be unity. 
The characteristics of the mixed finite element mod- 
els used herein are summarized in table 1. No conver- 
gence or other difficulties were encountered in apply- 
ing the foregoing computational procedure to the two 
problems. The converged solutions obtained by the 
mixed models were found to be in close agreement 
with the analytic solutions of Updike and Kalnins 
(1972) and the finite element solutions of Stein, 
Wagner, and Wriggers (1986) and Simo et al. (1984). 
Typical results are presented in figures 4 through 10 
for the hemispherical shell and in figures 12 through 
15 for the circular ring. The results presented in fig- 
ures 4 through 15 are considerably more extensive 
than previously reported and are intended to give 
physical insight into the response of the two struc- 
tures. All the numerical studies were performed on 


the Control Data Corp. CYBER 180-860 computer 
at NASA Langley Research Center. 

Hemispherical Shell in Contact With a Rigid 
Plate 

The first problem considered is the hemispherical 
shell pressed against a rigid flat plate. The geometric 
and material characteristics of the shell are shown 
in figure 3. Two cases are considered. In the first 
case the shell is subjected to a uniform ring load P 
and in the second case the shell is subjected to 
uniform internal pressure po in addition to the ring 
load. The first case is similar to that presented by 
Stein, Wagner, and Wriggers (1986) and Updike and 
Kalnins (1972). However, even in this case the results 
presented herein are for a much higher loading range 
than that considered in the cited references. 

The shell was modeled using uniform and non- 
uniform grids of the M3-2 and M4-3 elements (see ta- 
ble 1), with small-size elements in the contact zone. 
The converged solutions presented herein correspond 
to a nonuniform grid of 30 M3-2 elements. The char- 
acteristics of the nonuniform grids used are sum- 
marized in table 2. The variations in the vertical 
displacement at the loaded end and in the strain en- 
ergy components with the ring loading are shown in 
figure 4. The deformed configurations of the shell 
for different combinations of P and po are shown in 
figure 5. The distribution of the strain energy densi- 
ties and the contact pressures along the meridian, for 
different values of po and P, are shown in figures 6 
and 7. As can be seen from figure 6, both the bending 
and the transverse shear strain energy has a localized 
character. The variation of the contact region with 
loading is shown in figure 8. In the absence of in- 
ternal pressure a gap develops between the shell and 
the plate, with increasing loading. However, the size 
of the gap is too small to be noticeable (see fig. 5), 
and, as expected, the presence of internal pressure 
reduces the gap and increases the contact area. For 
Po/E = 5 X 10”^, no gap develops throughout the 
range of the ring loading considered. The predictions 
of the M3-2 and M4-3 models are shown in figure 9. 
For the same number of degrees of freedom, the two 
models result in comparable accuracy. This may be 
attributed to the fact that the accuracy is governed 
by the contact conditions. 

The number of iterations required for both the 
load and the arc-length incrementation is given in 
table 3. The number of iterations required for the 
two-level iteration scheme proposed by Wriggers and 
Nour-Omid (1984) is also given in table 3. As 
can be seen from table 3, use of the arc-length 
incrementation in conjunction with the single-level 
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iteration scheme is more efficient than the other 
schemes. 

The effect of the magnitude of the penalty pa- 
rameter on the accuracy of the total strain energy 
obtained by the mixed models, at Pto/Eh^ = 2.0 
and po = 0, is depicted in figure 10. As can be seen, 
the accuracy of the total strain energy obtained by 
the mixed models is fairly insensitive to the choice of 
e in the range of efEh from 10^ to 10^®. 

Circular Ring Pressed Against a Rigid Plate 

The second problem is a circular ring pressed 
against a rigid plate by a concentrated load at the 
top. The geometric and material characteristics of 
the ring are shown in figure 11. The ring was modeled 
by using uniform grids of M3-2 and M4-3 elements. 

The deformation of the ring is nearly in- 
extensional, with both the extensional and the trans- 
verse shear strain energy equal to zero throughout 
the range of loading considered. The total strain en- 
ergy is essentially equal to the bending energy. The 
accuracy of the vertical displacement Wa and the 
strain energies U obtained by using the M3-2 and 
M4-3 models is depicted in figure 12. Again, both 
models gave comparable accuracy. The deformed 
configurations of the ring are shown in figure 13 along 
with the distribution of the total strain energy den- 
sity along the meridian, for different values of the 
concentrated load. 

The variation of the contact region with loading is 
depicted in figure 14. As can be seen from figure 14, 
only a small segment of the ring remains in contact, 
for any given value of the load. However, the gap is 
very small and cannot be noticed in figure 13. 

The effect of the magnitude of the penalty pa- 
rameter on the accuracy of the total strain energy 
obtained by the mixed models, at Pr^jEI = 2.863, 
is depicted in figure 15. Again the accuracy of the 
total strain energy is fairly insensitive to the choice 
of e in the range of e/EA from 10^^ to 10^. 

Concluding Remarks 

Simple mixed finite element models and a compu- 
tational procedure are presented for the solution of 


frictionless contact problems. The analytical formu- 
lation is based on a form of Reissner’s large-rotation 
theory of the structure with the effect of transverse 
shear deformation included. The contact conditions 
are incorporated into the formulation by using a per- 
turbed Lagrangian approach with the fundamental 
unknowns consisting of the internal forces (or stress 
resultants), the generalized displacements, and the 
Lagrange multipliers associated with the contact con- 
ditions. The elemental arrays are obtained by using 
a modified two-field mixed variational principle. The 
modification consists of augmenting the functional of 
that principle by two terms: the Lagrange multiplier 
vector associated with the nodal contact pressures 
and a regularization term which is quadratic in the 
Lagrange multiplier vector. 

The shape functions used in approximating the 
generalized displacements and the Lagrange multi- 
pliers are selected to be the same and are different 
from those used in approximating the internal forces. 
The internal forces and the Lagrange multipliers are 
allowed to be discontinuous at interelement bound- 
aries. The nonlinearities due to both the large rota- 
tions and the contact conditions are combined into 
the same iteration loop and are handled by using the 
Newton-Raphson iterative scheme. 

Two numerical examples of axisymmetric defor- 
mations of a hemispherical shell and planar deforma- 
tions of a circular ring are presented. Both structures 
are pressed against a rigid plate. The detailed infor- 
mation presented herein is considerably more exten- 
sive than previously reported and helps in gaining 
physical insight about the response of both struc- 
tures. The numerical studies have demonstrated the 
high accuracy of the mixed models and the effective- 
ness of the computational procedure, which combines 
both the geometrically nonlinear terms and the con- 
tact conditions in one iteration loop. 

NASA Langley Research Center 
Hampton, VA 23665-5225 
February 6, 1989 
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Appendix A 

Fundamental Equations of the 
Large-Rotation Theories Used 
in the Present Study 

The fundamental equations of the large-rotation 
theories of Reissner (1972a, 1972b) for shells of rev- 
olution and curved beams are summarized herein. 
For shells of revolution, only axisymmetric deforma- 
tions are considered, and the effects of transverse 
shear deformation and moments turning around the 
normal to the middle surface are included. For 
curved beams, only planar deformation is considered, 
and the effects of transverse shear deformation and 
centerline extensibility are included. A total La- 
grangian description of the deformation is used, and 
the structure configurations at different load levels 
are referred to the initial coordinate system of the 
undeformed structure. 

Strain-Displacement Relationships 

The relationships between the strains and the 
generalized displacements are given as follows. For 
both the shell of revolution and the curved beam, 
the expressions for £5,^5, and 75 are 

£3 = cos (0 + (f)o) dgU — sin(</) + <l>o) dgW + cos (j) — I 

= (1 + £50) cos 7o - 1 

Ks = ds(f> 

75 = sin((/) -h (/>o) + cos(0 -f (f>o) dsw -h sin</> 

= (1 -h£5o)sin7o 


In addition, for the shell of revolution, the expres- 
sions for £^, and K>n are 

u 


- [sin(0 -h <t)o) - sin 0o] 
r 

«n = - [cos(<^ + <l>o) - COS 0o] 
r 

where Sgo is the extensional strain (relative change of 
length in the 5-direction); 70 is the transverse shear- 
ing strain; Ks is the bending strain in the 5-direction; 
6s and 75 are virtual extensional and shearing strains; 
E0 and kq are the extensional and bending strains in 
the circumferential directions for the shell; Kn is the 
bending strain associated with the moment turning 
about the normal to the middle surface of the shell; 
and ds = dfds. 

Constitutive Relations 

The material of the structures is assumed to be 
linearly elastic and isotropic. The relations between 
the strain components and the internal forces (stress 
resultants) are given in table Al. In table Al, 
and v are Young’s modulus, shear modulus, 
and Poisson’s ratio of the material; h is the shell 
thickness; and A and / are the cross-sectional area 
and the moment of inertia of the curved beam. 


Table Al. Constitutive Relations 



Shell of revolution 

Curved beam 

Extensional group 


f 1 _ 1 r 1 -v' 

\ eg J Eh —1/ 1 

(Ns' 

{Ng^ 



Bending group 


12 

' \ -V 

1 


fM.| 

M. 

1 


~ Eh^ 

-V 1 


< 



1 . 


[ • • 1- 



[m„ J 

Transverse shear 







Is = ^Qs 
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Appendix B 

Formulas for the Elemental Arrays [Q], [/2], 
and [go] 

The explicit forms of the elemental arrays [Q], [/?], 
and {qo} are given in this appendix. For convenience, 
each array is partitioned into blocks according to 
contributions from displacement and contact nodes. 
The expressions of the typical partitions (or blocks) 
are given in table Bl. In table Bl, AT* and 
are the shape functions for the Lagrange multipliers 


and generalized displacements; m is the number of 
displacement nodes in the element; c is the number 
of nodal points in contact in the element; and is 
the element domain. The range of the indices i and 
j is from 1 to c, and the range of the index i' is from 
1 to m; and <g>is the unit ramp (or singularity) 
function defined as follows: 


< 9 > 




r 


where g ^ —g and n = 0 or 1. 


( 5 > 0 ) 

(9 < 0 ) 


Table Bl. Explicit Form of Typical Partitions of 
Arrays [Q], [R], and {^o} 


Array 

Number of 
partitions or blocks 

Formula for typical partition 

[Q] 

m X c 

n^'n^ <g>° an 

[R] 

cx c 

- <g>^ dn 

{9o} 

c 

Jfl(e) ^ ^ go ^ d,Cl 
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Appendix C 

Computational Procedure Based on Generalized Arc-Length Control 

The computational procedure (see Ramm 1981) used in the present study is summarized in this appendix. 

Preprocessing and Initial Calculation Phase 

• Step 1. Input problem and model data; select values for the parameters a and select estimates for the 
penalty parameter and the generalized arc- length increment AC (see Noor and Peters 1981a, 1981b); and 
assume the contact status at the different nodes. 

• Step 2. Generate linear elemental arrays. 


Solution Phase 

• Begin incrementation loop. 

• Begin Newton-Raphson iteration loop. 

• Step 3. Predict the generalized displacements, internal forces (stress resultants), and Lagrange 
multipliers for the structure. 

• Step 4. Generate the vector {dZ/dQi using the following set of equations; 



_ 

r^l 

dp 

IflC/i 

1 dp J 

U dC. 


where the range of I and J is from 1 to the number of degrees of freedom in the model; and subscript 
i refers to the rth incrementation step. The positive sign for dpfdC is taken for the stable equilibrium 
path and the negative sign for the unstable path (see Noor and Peters 1981a, 1981b). 

• Step 5. Complete initial estimates for (Ap)j and {AZ}j as follows: 


(Ap)f> = |[AC 

{AzW)j = {^}^ AC 

• Step 6. Generate nonlinear elemental array; eliminate the internal forces and the Lagrange multipliers 
from the elemental equations; and assemble the left- and right-hand sides of the equations. 

• Step 7. The change {AZ} in the response vector {Z} during the ith incrementation step can be 
expressed as the sum of the two vectors 


{AZ}i = {AZ}i + (Ap)i{AZ}i 


The vectors {AZ} and {AZ} are obtained by solving the following system of algebraic equations: 


[K] + 


W 
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where superscript r refers to the rth iteration cycle, and subscript i has been dropped for convenience. 
I • Step 8. Calculate the incremental load factor Ap using the equation 


Ap^*") = 


-a{AZW}^{A^’'^} 

a {aZ(i)}‘ {azW} +/3Ap(l) 


• Step 9. Update the load factor and the response vector 


j p(’-+l) = p(’-) + Ap(r) 

= {zW} + {azW} 

• Step 10. Check the contact status (and modify the contact conditions) at each node as needed: If 

; p > 0 and A > 0, then the constraint is inactive; if p < 0, then the constraint is active. If contact 

( status is the same as previously assumed, then continue. Otherwise, add the contact contributions 

for nodes with active constraint and delete these contributions for nodes with inactive constraint; 
then return to step 6. 

i 

I • Step 11. Check convergence of Newton- Raphson iterations, that is, 

n 

^ where n is the total number of degrees of freedom in the model; and the tolerance is prescribed. If 

I convergence is achieved then continue; otherwise return to step 6. 

f 

j • Step 12. If p > Pmaxj then stop; otherwise set C = C + and return to step 3. 


I 


I 
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Table 1. Characteristics of Mixed Finite Element 
Models Used in Numerical Studies 



Number of 

Maximum number 


Number of 


displacement 

of Lagrange 

Number of parameters 

quadrature 

Designation 

nodes 

multiplier nodes 

per internal force 

points® 

M3-2 

3 

3 

2 

3 

M4-3 

4 

4 

3 

4 


“All elemental arrays are evaluated using Gauss-Legendre quadrature formulas except for [Q], [fl], 
and {ffo} which are evaluated using Newton-Cotes formulas. 


Table 2. Characteristics of Nonuniform Grids Used in Modeling Hemispherical Shell 

[Element length I = 


8 elements 

12 elements 

15 elements 

30 elements 

Element 


Element 

^<Po, 

Element 

^4>o, 

Element 


no. 

deg 

no. 

deg 

no. 

deg 

no. 

deg 

1 

1.5 

1 

1.0 

1 to 10 

3 

1 to 20 

1.5 

2 to 4 

9.5 

2 to 6 

5.8 

11 

5 

21 

2.0 

5 to 8 

15.0 

7 to 12 

10.0 

12 

9 

22, 23 

5.0 





13 

13 

24 

6.0 





14 

16 

25 to 30 

7.0 



1 


15 

17 
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I 


I 


Table 3. Comparison of Number of Iterations Required for Present Single-Level Iteration 
Procedure With Two-Level Iteration Procedure 


Hemispherical shell in contact with a rigid plate (see fig. 3); 
15 M3-2 elements along the meridian 



Arc-length incrementation 

Load incrementation 



Single 


Two 


Single 

Two 


Pro 

Eh'^ 

iteration 

Pro 

Eh^ 

iteration 

Pro 

Eh^ 

iteration 

iteration 

Step 

loop® 

loops^ 

loop^ 

loops^ 

1 

2.3442 

5(3) 

2.3442 

7 

4.0 

8(5) 

17 

2 

9.6426 

5(2) 

9.6411 

8 

8.0 

7(6) 

21 

3 

15.315 

5(10) 

15.318 

17 

12.0 

8(5) 

21 

4 

21.927 

5 (5) 

21.939 

14 

16.0 

6(3) 

17 

5 

29.316 

7(5) 

29.349 

14 

20.0 

6(7) 

8 

6 

36.775 

7(4) 

36.845 

9 

24.0 

4(7) 

11 

7 

44.466 

6(4) 

44.541 

9 

28.0 

4(3) 

10 

8 

51.627 

4(8) 

51.818 

5 

32.0 

3(8) 

11 

9 





36.0 

4(3) 

5 

10 





40.0 

4(5) 

5 

11 





44.0 

8(8) 

21 

12 





48.0 

3(3) 

7 

13 


■■ 



52.0 

7(3) 

18 

Total .... 

44 (41) 


83 

i 


72 (66) 

172 


“Numbers in parentheses indicate the number of iterations for the uniform grid; other numbers are for 
the nonuniform grid (see table 2). The prescribed tolerance was selected to be 1 x 10“^. 

*Based on the solution procedure of Wriggers and Nour-Omid (1984). 
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Figure 1. Curved beam and shell elements and sign convention. 


Rigid plate 




Tn 



I Contact region 
: and contact 

1 pressure 
! distribution 


Figure 2. Characteristics of frictionless contact of a structure pressed against a rigid plate. 



Rigid plate 



E = 21 X 10"'®Pa 
V =0.5 
rQ = 1.0 m 

h = 3.33x 10‘^m 


Boundary conditions 
At b: (|) = 0 


Figure 3. Hemispherical shell in contact with a rigid plate. 




Figure 4. Variation of vertical displacement at the loaded end and total strain energy with loading for 
hemispherical shell. 







Figure 7. Distribution of contact pressures along the meridian of hemispherical shell. 



s/h s/h 

Figure 8. Variation of the contact region with loading for hemispherical shell. Sj and So refer to s-coordinates 
of initial and last contact points. 




El = 0.1 N-m'^ 



EA = 10^ N 
GA= 10^ N 
ro= 1.0 m 


Figure 11. Circular ring pressed against rigid plate. 


Finite 




Figure 12. Variation of vertical displacement at the load point and total strain energy with loading for circular 
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Prjf/EI Wg/ro 
1.464 - 0.6 




Figure 13. Deformed shapes and total strain energy density distributions for circular ring. 



Figure 14. Variation of contact region with loading for circular ring. Sj and Se refer to s-coordinates of initial 
and last contact points. 
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